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Abstract 

The effect of reducing the formal order of accuracy of a finite-difference scheme in order to 
optimize its high-frequency performance is investigated using the 1-D nonlinear unsteady inviscid 
Burgers’ equation. It is found that the benefits of optimization do carry over into nonlinear appli- 
cations. Both explicit and compact schemes are compared to Tam and Webb’s explicit 7-point 
Dispersion Relation Preserving scheme as well as a Spectral-like compact scheme derived follow- 
ing Lele’s work. Results are given for the absolute and L2 errors as a function of time. 

Introduction 

Computational aeroacoustics is concerned with the time-accurate solution of flow and acous- 
tic phenomena over long periods of time. To accomplish this goal, high-order finite-difference 
schemes and optimized schemes have been developed [e.g.. Refs 1-8]. 

There are two main classes of high-accuracy finite-difference schemes: explicit schemes and 
compact schemes. Explicit schemes employ large computational stencils for accuracy, while com- 
pact schemes use smaller stencils by using the flux derivatives as independent variables at each 
grid point. While compact schemes are more accurate than the equivalent explicit scheme, solving 
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for each flux derivative requires a scalar tridiagonal or pentadiagonal matrix inversion 1 . Recently, 
however, an new prefactorization method has been introduced that allows the matrices to be pref- 
actored into two easier to solve matrices 7,8 . 

There are also two strategies for designing finite-difference stencils. For a given stencil, one 
can either choose coefficients that return the highest formal order of accuracy possible, or the 
stencil can be optimized to reduce high-frequency errors at the price of lowering the formal order 

of accuracy. This optimization work has been used for both compact 1 and explicit 3 stencils. 

In this work, the nonlinear performance of optimized schemes will be investigated using the 1 - 
D unsteady Burgers’ equation. 

Mathematical and Numerical Formulation 

Given a function f(x), a general finite-difference central derivative at grid point n can be writ- 
ten as: 
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The right-hand-side of Eq. (1) determines if the differencing scheme is explicit or compact. 
For a compact differencing scheme, each derivative depends on the value of its neighboring deriv- 
atives, thus, a scalar matrix inversion is required to obtain the values of the derivatives. For an 
explicit differencing scheme, however, the value of each derivative is independent of its neigh- 
bors. 

The error term is a function of the order of accuracy of the scheme. To estimate the magnitude 
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of the error, we follow the work of Lele 1 and Tam and Webb 3 . Let us assume that the function we 
are taking the numerical derivative of is a sine wave of period one: 


/ = sin(27tx) (2) 

— = 2kcos(2kx) 


where -0.5 < x < 0.5. We have a given number of points per wavelength: 


Ax = 


1 

PPW 


(3) 



Figure 1 : Error Magnitude as a Function of Points Per Wavelength for Finite-Difference Deriva- 
tives 
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Thus, the function at the grid points can be written as: 


/,• = sin(2ju'Aj:) 


(4) 


and the numerical derivative at x = 0 becomes: 
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The coefficients for the schemes investigated are given in Table 1. 

The error in a central derivative is dispersive, and the numerical derivative takes the form: 



Figure 2: Order of Accuracy as a Function of Points Per Wavelength for Finite-Difference Deriv- 
atives 
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Using this definition, the derivative at x=0 can be written as: 
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The error at x = 0 can then be written as: 
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Table 1: Finite-Difference Scheme Coefficients 


Scheme 

bi 

a l 

a 2 

a 3 

Order of 
accuracy 

explicit 4 

0 

2/3 

1/12 

0 

4 

7-point 

DRP 

0 

0.7708824 

-0.1667059 

0.0208431 

4 

explicit 6 

0 

3/4 

-3/20 

1/60 

6 

compact 4 

1/4 

3/4 

0 

0 

4 

optimized 

compact 

0.35196 

0.7839867 

0.0339867 

0 

4 

compact 6 

1/3 

7/9 

1/36 

0 

6 
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Figure 1 shows the error magnitude as a function of the number of points per wavelength. 
Notice the wide range of errors for different schemes even though they have the same order of 
accuracy. Figure 2 shows the order of accuracy as a function of the number of points per wave- 
length. Notice that it requires 25-40 points per wavelength for the schemes to settle to their for- 
mal order of accuracy. In practice, this many points cannot be used due to computer limitations. 

In Fig. 1, the two optimized schemes show similar error curves. The effect of the optimization 
is to increase the resolution of the scheme with few points per wavelength by lowering the formal 
order of accuracy, and hence reducing the performance of the scheme with many points per wave- 
length, The curves in Fig. 1 show that the optimized schemes always perform better than the 
unoptimized scheme of the same order of accuracy (while also requiring more work). The curves 
also show that there is a range where the optimized scheme will be more accurate than a higher- 
order scheme using the same stencil, but this range is fairly limited. 

For a linear problem with a small range of frequencies, the optimized scheme can be made to 
outperform the higher-order scheme. However, with a change in grid spacing, the higher-order 
scheme can also be made to outperform the optimized scheme. As the frequency range increases, 
it becomes less clear which scheme provides a better solution. 

In an umsteady nonlinear problem, however, the frequency range changes as time progresses. 
An example of a nonlinear equation is the inviscid 1-D Burgers’ equation: 


du d_ 
dt + dx 



= 0 


(9) 
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Given an initial single-frequency periodic solution on the domain 0 < x < 1: 


u(x,0) = c 0 + Asin(27tjc), 


( 10 ) 


the solution can be written as: 


u(x,t) = k( 5,0) = c 0 + Asin(2it^) 


( 11 ) 


where 


5 = x-ut (12) 

0 <£<! 


However, at time: 


t = 


1 

2 7t A 


(13) 


a shock forms as the wave steepens, and the shock is located using Whitham’s area rule 9 : 


27tA sin0 


(14) 
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(15) 


(16) 

= c 0 + Asin(2^j) 

(17) 

= c 0 + Asin(27t^ 2 ) 

(18) 


The remaining part of the solution is found using Eq. (11). 

Figure 3 shows a solution of Burgers’ equation, with the initial condition: 

u(x,0) = | + isin(2TCx) (19) 



Figure 3: Solution of 1-D Nonlinear Burgers’ Equation With Single-Frequency Input 



In this solution, the wave steepens as it travels, forming a shock at time = 0.637. As time pro- 
ceeds, the lower-frequency waves transfer more and more energy to the high-frequency waves, as 
illustrated in Fig. 4. 

Since this is the inviscid Burgers’ equation, the shock has no thickness. This means that there 
is no upper limit to the frequencies in the solution once the shock appears. Numerically, this 
means that there is no grid spacing that will resolve the shock. 

Since this equation allows unlimited frequency growth and has an exact solution, it was used 
to compare the nonlinear performance of the six schemes given in Table 1. 

Numerical Test Description 

Two numerical tests will be shown. Test 1 has the initial distribution: 



Figure 4: Frequency Content as a Function of Time for Solution of Burgers’ Equation 
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( 20 ) 


u(x , 0) = \ + ^sin(2j:x) 
4 4 


while Test 2 has the initial distribution: 


u(x,0) = | + ^sin(27tx) (21) 


In Test 2, the lower amplitude of the sine wave reduces the rate of growth of the high-fre- 
quency modes, causing the wave to steepen more slowly. This lesser rate of nonlinearity will 
allow more time for any advantages of optimization to be shown before the high-frequency waves 
move into the unresolved range. 

The grid for this problem has 24 equally spaced points, giving 24 points per wavelength for 
the initial wave. From Fig. 2, this corresponds roughly to the minimum number of points for each 
scheme to obtain its formal order of accuracy. This is also more than the number of points usually 
used for linear wave propagation in computational aeroacoustics calculations. The domain is peri- 
odic in order to remove the effect of boundary stencils and boundary condition specification from 
the results. The periodic domain may favor the compact schemes more than the explicit schemes, 
due to the effects of boundary stencil specifications on compact schemes 10 . 

The time stepping method used is Stanescu and Habashi’s 11 fourth-order nonlinear extension 

of Hu’s 5-6 optimized Runge-Kutta method 12 . The time step used is CFL = 0.025. The small time 
step is used to ensure that the error is from the spatial discretization and not the time marching 
scheme. 
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The compact differences are prefactored using Hixon and Turkel’s method 7 , and calculated 
using the method described by Hixon 8 . The forward and backward stencils for the optimized com- 


pact scheme are given as: 

( 1 - a)Z>f + aDf + 1 = P/, + , + ( 1 - 2p)/, + O - 1)/,_ , (22) 

(1 -a)Df + aZ>f_, = (l-p)/ j + 1 + (2P-l)/ 1 .-p/ 1 .. 1 


where 


a = 0.29157494 (23) 

(3 = 0.86318328 

Results 

Calculations were performed for Tests 1 and 2. No filtering or artificial dissipation terms were 
added; while, the solutions would be improved by filtering out unresolved frequencies, it was 
found that the error due to the filtering overwhelmed the errors from the schemes as the shock 
formed. The reason was that the 10th and 12th order explicit filters that were tested damped high 
frequency components of the solution that the more accurate compact schemes could resolve. 

Results are shown for the error amplitude as a function of time for each numerical scheme. 
Both the L 2 and max error norms are shown. It was found that the max norm favored the opti- 
mized schemes, while the L 2 norm was weighted more toward the high-order schemes. Since the 
L 2 norm contains both the low-frequency and high-frequency error components, it was considered 
to be the more accurate measure of error for this problem. 
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Figure 5: Time History of the L2 Norm Error Amplitude for Test Input 1 
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Figure 6: Time History of the Max Norm Error Amplitude for Test Input 1 
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Figure 7: Time History of the L2 Norm Error Amplitude for Test Input 2 



time 

Figure 8: Time History of the Max Norm Error Amplitude for Test Input 2 
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Figures 5 and 6 show the error time history for Test 1. In both figures, we see the initial advan- 
tages of the high-order schemes. However, at time = 0.2, the wave has steepened to a point that 
high-frequency components of the solution exist. After this point, the optimized schemes perform 
equal to or slightly better than the high-order schemes with the same stencil size. This is more 
apparent using the max norm of error, but also is shown in the L 2 norm of the error. 

Figures 7 and 8 show the error time history for Test 2. Notice that the shock takes much longer 
to form in this test. In Figs. 7 and 8, the same behavior is shown as in Test 1, with the optimized 
schemes showing more error reduction. This is due to the slower growth of high-frequency com- 
ponents in this test. Figure 8 shows that the max error norm is significantly lower with the opti- 
mized schemes, while Fig. 7 shows that the improvement is still seen in the L 2 norm results. 

Conclusions 

A nonlinear comparison of optimized and high-order finite-difference schemes has been per- 
formed using the 1-D inviscid Burgers’ equation. From these simple tests, several conclusions can 
be drawn. 

First, the compact schemes perform better than the explicit schemes, usually giving an order 
of magnitude reduction in the errors. 

Second, the formal order of accuracy of a scheme is much less important than its performance 
at marginal resolution. Again, this is shown by the performance of the compact schemes and opti- 
mized schemes for these test problems. 

Third, in most realistic nonlinear calculations, there will be areas of marginal resolution in the 
grid. The error due to these areas will dominate the error in the solution. The improved perfor- 
mance of the optimized schemes at marginal resolution can more than compensate for the slightly 
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increased error in fully resolved regions. 

In summary, while these results are necessarily not general, they do indicate that the opti- 
mized schemes show an advantage for nonlinear problems as well as linear problems due to their 
improved performance at marginal resolution. Since the high-frequency waves are the major 
source of error, the gain from optimization more than offsets the loss of formal accuracy for both 
explicit and compact differencing schemes in these tests. 
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